% Script to predict local magnetic declination variation from 
% local F, and remote XYZ at two obs


%% load data.

basedir = '/Users/manojnair/data/obs_mag_data/Canada/Absolute_netCDF/';
%basedir = '/Users/manojnair/data/obs_mag_data/Canada/Secular_Removed_netCDF/';

start_index = 0; % start from the beginning
ndata = 9468000; % read all data

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));

ncfname = [basedir 'TAL_1995_2012_Absolute.nc'];

[x_TAL, y_TAL, z_TAL , X_ID, Y_ID, Z_ID, obj_TAL] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

ncfname = [basedir 'GHC_1995_2012_Absolute.nc'];
[x_GHC, y_GHC, z_GHC , X_ID, Y_ID, Z_ID, obj_GHC] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

ncfname = [basedir 'PBC_1995_2012_Absolute.nc'];
[x_PBC, y_PBC, z_PBC , X_ID, Y_ID, Z_ID, obj_PBC] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

ncfname = [basedir 'CBB_1995_2012_Absolute.nc'];
[x_CBB, y_CBB, z_CBB , X_ID, Y_ID, Z_ID, obj_CBB] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

% Interpolate data points with gaps less than N. This should be run only once. 
maxgap = 11; 

x_TAL = shortinterpolate(x_TAL,maxgap);
y_TAL = shortinterpolate(y_TAL,maxgap);
z_TAL = shortinterpolate(z_TAL,maxgap);

x_GHC = shortinterpolate(x_GHC,maxgap);
y_GHC = shortinterpolate(y_GHC,maxgap);
z_GHC = shortinterpolate(z_GHC,maxgap);

x_PBC = shortinterpolate(x_PBC,maxgap);
y_PBC = shortinterpolate(y_PBC,maxgap);
z_PBC = shortinterpolate(z_PBC,maxgap);

x_CBB = shortinterpolate(x_CBB,maxgap);
y_CBB = shortinterpolate(y_CBB,maxgap);
z_CBB = shortinterpolate(z_CBB,maxgap);


%% find the common data with GHC and TAL


L = isnan(x_TAL) | isnan(y_TAL) | isnan(z_TAL) | ...
    isnan(x_GHC) | isnan(y_GHC) | isnan(z_GHC) | ...
    isnan(x_PBC) | isnan(y_PBC) | isnan(z_PBC) | ...
    isnan(x_CBB) | isnan(y_CBB) | isnan(z_CBB);

L = ~L;
%% prepare input data with long continous segments
% idea from http://stackoverflow.com/questions/2212201/


D = diff(L);

b.beg = 1 + find(D == 1);

if L(1)
  b.beg = [1;b.beg];
end

b.end = find(D == -1);

if L(end)
  b.end(end+1) = numel(L);
end

seg_length = b.end - b.beg;

seg_length_sorted = sort(seg_length);
%% get the longest continuous seg

k = find(seg_length == 17366);

[input, output] = prepareanndata(k, b, time_array_min, x_TAL, y_TAL, z_TAL, x_GHC, y_GHC, z_GHC, x_PBC, y_PBC, z_PBC, x_CBB, y_CBB, z_CBB, ...
  obj_TAL, obj_GHC, obj_PBC, obj_CBB);

k = find(seg_length == 11973);

[input1, output1] = prepareanndata(k, b, time_array_min, x_TAL, y_TAL, z_TAL, x_GHC, y_GHC, z_GHC, x_PBC, y_PBC, z_PBC, x_CBB, y_CBB, z_CBB, ...
  obj_TAL, obj_GHC, obj_PBC, obj_CBB);

%% plot data

subplot(211);

plot(output,'r');
hold;
plot(input(:,1),'b');
title(sprintf('rms = %5.2f', rms(output-input(:,1))));
subplot(212);
plot(output1,'r');
hold;
plot(input1(:,1),'b');
title(sprintf('rms = %5.2f', rms(output1-input1(:,1))));


